In-class 8 Nov
Session preparation
# Setup ---------------------...
library(leaflet)
library(tmap)
library(sf)
library(tidyverse)
# Get the data:
dc <-
st_read('data/raw/rasters/dc.shp')## Reading layer `dc' from data source
## `/Users/brianevans/Library/Mobile Documents/com~apple~CloudDocs/ppol/data/raw/rasters/dc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -77.11976 ymin: 38.79164 xmax: -76.90939 ymax: 38.99584
## Geodetic CRS: NAD83dc_crimes <-
read_csv('data/raw/dc_crimes.csv')
census <-
st_read('data/raw/census/census.shp')## Reading layer `census' from data source
## `/Users/brianevans/Library/Mobile Documents/com~apple~CloudDocs/ppol/data/raw/census/census.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3492 features and 17 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -83.67539 ymin: 36.54085 xmax: -74.98628 ymax: 39.72304
## Geodetic CRS: NAD83rasters <-
read_rds('data/raw/rasters/dc_lc.rds')
# Some data wrangling -------...
# Subset dc_crimes to violent crimes:
violent_crimes <-
dc_crimes %>%
filter(lubridate::year(date_time) == 2020,
offense_group == 'violent')
# Convert violent_crimes to an sf object:
crimes_sp <-
st_as_sf(
violent_crimes,
coords = c('longitude', 'latitude'),
crs = 4326)
# Subset census data to DC and columns of interest:
dc_census <-
census %>%
filter(STATEFP == 11) %>%
select(GEOID, edu, income, population)
# Convert census and crimes data to EPSG 5070:
list(crimes_sp,
dc_census) %>%
map(~st_transform(., crs = 5070)) %>%
set_names('crimes_prj', 'census_prj') %>%
list2env(.GlobalEnv)## <environment: R_GlobalEnv># Join crimes and census data, then simplify the shapes:
census_simple <-
st_join(
census_prj,
crimes_prj) %>%
rmapshaper::ms_simplify(keep = 0.05)
# Calculate the number of crimes per 1000 residents and
# return as a simple features object:
census_simpler <-
census_prj %>%
left_join(
census_simple %>%
as_tibble() %>%
group_by(GEOID) %>%
summarize(crimes_per_1000 = n()/unique(population)*1000),
by = 'GEOID')
# Transform shapes to 4326:
list(crimes_prj, census_simpler) %>%
map(~ st_transform(., crs = 4326)) %>%
set_names('crimes_4326', 'census_4326') %>%
list2env(.GlobalEnv)## <environment: R_GlobalEnv># Subset crimes to robberies:
robberies <-
crimes_4326 %>%
filter(offense == 'robbery')
# Rasterize crimes:
crimes_raster <-
crimes_prj %>%
dplyr::select(-everything()) %>%
st_transform(crs = st_crs(rasters)) %>%
as_Spatial() %>%
raster::rasterize(
rasters, # These are the cells that the points are rasterized to
fun = 'count',
background = 0)Homework
Question One: Raster reclassification
Open water is classified as the value 11 in the NLCD raster. Using a logical statement, classify the raster such that all open water raster cells are given the value 0 and all non-water cells are given the value 1. Assign the raster to your global environment with the name land.
land <-
rasters$nlcd != 11Modify your land raster such that all water raster cells are given the value NA and all non-water raster cells are given the value 1.
land[land == 0] <-
NAPlot the resultant land raster:
raster::plot(land)
Question Two: Urban landscapes
Create a reclass matrix to classify the impervious_surface raster into low, medium and high-intensity urban land cover classes.
tribble(
~ from, ~ to, ~becomes,
0, 10, 1,
10, 60, 1,
60, 100, 1) %>%
as.matrix()## from to becomes
## [1,] 0 10 1
## [2,] 10 60 1
## [3,] 60 100 1Use cut() to ensure that your reclass matrix will classify your raster as expected.
tibble(values = 1:10) %>%
mutate(
class =
cut(
values,
breaks = c(0, 1, 6, 10),
labels = c('low', 'medium', 'high'),
include.lowest = TRUE,
right = TRUE))## # A tibble: 10 × 2
## values class
## <int> <fct>
## 1 1 low
## 2 2 medium
## 3 3 medium
## 4 4 medium
## 5 5 medium
## 6 6 medium
## 7 7 high
## 8 8 high
## 9 9 high
## 10 10 highReclassify the impervious_surface raster based on your reclass matrix and assign the resultant object to your global environment with the name urban_intensity.
tribble(
~ from, ~ to, ~becomes,
0, 10, 1,
10, 60, 2,
60, 100, 3) %>%
as.matrix() %>%
raster::reclassify(
rasters$impervious_surface,
.,
include.lowest = TRUE,
right = TRUE)## class : RasterLayer
## dimensions : 747, 621, 463887 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 1610775, 1629405, 1913775, 1936185 (xmin, xmax, ymin, ymax)
## crs : +proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
## source : memory
## names : impervious_surface
## values : 1, 3 (min, max)Question Three: Aggregate:
In a single piped statement (prompts provided as comments):
crimes_raster %>%
# Aggregate the crimes_raster to a cell size of roughly 500 x 500 m:
raster::aggregate(
fact = 17,
fun = sum) %>%
# Mask the data to the dc shapefile:
raster::mask(
st_transform(
dc,
crs = st_crs(crimes_raster))) %>%
# Plot the data using ggplot():
raster::rasterToPoints() %>%
as_tibble() %>%
rename(crimes = layer) %>%
ggplot() +
geom_tile(aes(x = x, y = y, fill = crimes)) +
scale_fill_viridis_c(option = 'viridis') +
coord_equal() +
theme_void()
Question Four
Generate an interactive map where each offense is given its own layer.
tmap_mode('view')
# Census tracts, colored by violent crimes per 1000 residents:
tm_shape(
name = 'census tracts',
census_4326) +
tm_polygons(
title = 'violent crimes per 1000 residents',
col = 'crimes_per_1000',
alpha = 0.5) +
# Homicides:
tm_shape(
name = 'homicide',
filter(crimes_4326,
offense == 'homicide')) +
tm_dots(size = 0.05,
clustering = TRUE) +
# Assault with dangerous weapon:
tm_shape(
name = 'assault w/dangerous weapon',
filter(crimes_4326,
offense == 'assault w/dangerous weapon')) +
tm_dots(size = 0.05,
col = 'blue',
clustering = TRUE) +
# Sex abuse:
tm_shape(
name = 'Sex abuse',
filter(crimes_4326,
offense == 'sex abuse')) +
tm_dots(size = 0.05,
col = 'red',
clustering = TRUE) +
# Robberies:
tm_shape(
name = 'Robbery',
filter(crimes_4326,
offense == 'robbery')) +
tm_dots(size = 0.05,
col = 'gray50',
clustering = TRUE)Leaflet
We’re going to generate an interactive plot using Leaflet.
We’re going to generate a map of the following files:
- crimes_4326
- census_4326
- An aggregated raster file
The basic map
We initiate our leaflet map (optimally) by choosing the file that best represents the spatial extent and applying the leaflet() function to the object:
census_4326 %>%
leaflet()From there, we can add layers. A good layer to add from the start is the basemaps, which here are called “tiles”:
census_4326 %>%
leaflet() %>%
addTiles()Adding and modifying polygons
We can use the addPolygons() function to add our shapes:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons()Our polygons are very ugly. Let’s make those lines thinner:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(weight = 2)And change the color:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor = 'yellow')To generate different colors representing value levels for the polygons, we need to specify a color palette function (see ?ColorBin):
population_pal <-
colorBin(
bins = 4,
palette = 'viridis',
domain = census_4326$population)And then define the palette as an argument inside of addPolygons():
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7)We can easily change the number of bins:
population_pal <-
colorBin(
bins = 8,
palette = 'viridis',
domain = census_4326$population)
# Add fill:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7)At this point, it’s useful to provide a legend. We specify the palette, the source for the values, and the opacity of the legend items:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7) %>%
addLegend(
pal = population_pal,
values = ~population,
opacity = 0.7)We can also change the position of the legend:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7) %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7)Adding points
We can add points as markers or circle markers. Here we’ll use markers, but see ?addMarkers.
To do so, if the map data are not defined using the point values, we need to extract the coordinates from the file (we’ll walk through this):
crimes_coords <-
st_coordinates(crimes_4326) %>%
as_tibble() %>%
bind_cols(as_tibble(crimes_4326))And then we specify the longitude and latitude columns inside of the addMarkers() function:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7) %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y) %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7)That’s an overwhelming amount of points! We can cluster the points using clusterOptions = markerClusterOptions():
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7) %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions()) %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 1)Now you!
Inside the addMarkers() function, use indexing to subset the markers to just offenses classified as “robbery”. Make the robberies circle markers (see ?addMarkers)
We can add labels to the map using the label = argument. Here, we’ll label the points with the date and time a violent crime occurred:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7) %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time) %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 1)Adding more basemaps
We can use the addProviderTiles() function to add basemaps from the leaflet community.
Type the following in your console to see a list of provider tiles:
providersHere, we’ll add digital orthophotos to our map:
census_4326 %>%
leaflet() %>%
addTiles() %>%
addProviderTiles(providers$Esri.WorldImagery) %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7) %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time) %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7)To see a list and demonstration of provider tiles go to this link: http://leaflet-extras.github.io/leaflet-providers/preview/index.html
Now you!
Use the link above to find and add a new provider tile to your map.
Layers control
You can provide users with the option to turn on or off layers. These come in two flavors:
- baseGroups: Controls and grouping of background tiles
- overlayGroups: Controls and grouping of all other layers
Let’s add a control for our baseGroups first. We define a group inside of our layer and use the addLayersControl() function to control the groups:
census_4326 %>%
leaflet() %>%
addTiles(group = 'Open Street Map') %>%
addProviderTiles(
providers$Esri.WorldImagery,
group = 'Orthophoto') %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7) %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time) %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7) %>%
addLayersControl(
baseGroups = c('Open Street Map',
'Orthophoto'))We do the same for points and polygons:
census_4326 %>%
leaflet() %>%
addTiles(group = 'Open Street Map') %>%
addProviderTiles(
providers$Esri.WorldImagery,
group = 'Orthophoto') %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7,
group = 'Population by Census Tract') %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time,
group = 'Violent Crimes') %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7) %>%
addLayersControl(
baseGroups = c('Open Street Map',
'Orthophoto'),
overlayGroups = c('Violent Crimes', 'Population by Census Tract'))We can also make an overlay group “born” hidden:
census_4326 %>%
leaflet() %>%
addTiles(group = 'Open Street Map') %>%
addProviderTiles(
providers$Esri.WorldImagery,
group = 'Orthophoto') %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7,
group = 'Population by Census Tract') %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time,
group = 'Violent Crimes') %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 1) %>%
addLayersControl(
baseGroups = c('Open Street Map',
'Orthophoto'),
overlayGroups = c('Violent Crimes', 'Population by Census Tract')) %>%
hideGroup('Violent Crimes')Now you!
Add your provider tile to the baseGroups and your robbery markers to the overlayGroups.
Raster layers
We can also add rasters to a Leaflet map. But first … we need to change the CRS of the raster to EPSG 4326 …
Now you!
In a single piped statement:
- Aggregate your
crimes_rasterto a resolution of roughly 500 x 500 m - Convert the raster to the same projection as
census_4326 - Mask the data to
census_4326 - Assign the object to your global environment with the name
raster_4326
We start by defining a palette for our raster object:
raster_pal <-
colorNumeric(
palette = 'Reds',
domain = raster::values(raster_4326))Then we add it with the function addRasterImage():
census_4326 %>%
leaflet() %>%
addTiles(group = 'Open Street Map') %>%
addProviderTiles(
providers$Esri.WorldImagery,
group = 'Orthophoto') %>%
addRasterImage(
raster_4326,
colors = raster_pal,
opacity = 0.9
) %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7,
group = 'Population by Census Tract') %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time,
group = 'Violent Crimes') %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7) %>%
addLayersControl(
baseGroups = c('Open Street Map',
'Orthophoto'),
overlayGroups = c('Violent Crimes', 'Population by Census Tract')) %>%
hideGroup('Violent Crimes')This is fairly awful looking. We can specify a color for NA inside of our palette:
raster_pal <-
colorNumeric(
palette = 'Reds',
domain = raster::values(raster_4326),
na.color = '#00000000')
census_4326 %>%
leaflet() %>%
addTiles(group = 'Open Street Map') %>%
addProviderTiles(
providers$Esri.WorldImagery,
group = 'Orthophoto') %>%
addRasterImage(
raster_4326,
colors = raster_pal,
opacity = 0.9
) %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7,
group = 'Population by Census Tract') %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time,
group = 'Violent Crimes') %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7) %>%
addLayersControl(
baseGroups = c('Open Street Map',
'Orthophoto'),
overlayGroups = c('Violent Crimes', 'Population by Census Tract')) %>%
hideGroup('Violent Crimes')We can add our raster to the overlayGroup as well:
raster_pal <-
colorNumeric(
palette = 'Reds',
domain = raster::values(raster_4326),
na.color = '#00000000')
census_4326 %>%
leaflet() %>%
addTiles(group = 'Open Street Map') %>%
addProviderTiles(
providers$Esri.WorldImagery,
group = 'Orthophoto') %>%
addRasterImage(
raster_4326,
colors = raster_pal,
opacity = 0.9,
group = 'Crime Density') %>%
addPolygons(
weight = 2,
fillColor =
~population_pal(population),
fillOpacity = 0.7,
group = 'Population by Census Tract') %>%
addMarkers(
lng = crimes_coords$X,
lat = crimes_coords$Y,
clusterOptions = markerClusterOptions(),
label = crimes_coords$date_time,
group = 'Violent Crimes') %>%
addLegend(
'bottomleft',
pal = population_pal,
values = ~population,
opacity = 0.7) %>%
addLayersControl(
baseGroups = c('Open Street Map',
'Orthophoto'),
overlayGroups = c('Violent Crimes',
'Population by Census Tract',
'Crime Density')) %>%
hideGroup('Violent Crimes')